import numpy as np
import random
from numpy import polyfit, poly1d


def designRs(T, Damp, ld, cd, fz, rarity):
    """
    # 地震波设计谱函数 用于求取一个频率在设计谱上对应的谱值
    # 根据《GB 50260 - 2013电力设施抗震设计规范》设计
    # 针对电力设施
    :param T: 特征周期序列
    :param Damp: 阻尼比
    :param ld: 地震烈度  6 7 7.5 8 8.25 8.5 9 10  对应到0-7
    :param cd: 场地类型 0 1 2 3 4
    :param fz: 场地分组 1 2 3
    :param rarity: 地震罕见度 frequent多遇地震 rare罕遇地震 basic基本地震动 very rare 极罕遇地震动
    :return:
    """
    # 根据抗震设防烈度，选择水平地震影响系数最大值
    arf_jg = 0
    # 水平地震影响系数最大值 6-0.05g 7-0.10g 7.5-0.15g 8-0.20g 8.25-0.25g 8.5-0.30g 9-0.40g 10-0.50g
    ACO_MAX = np.array([0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 1, 1.25])
    arf_max = ACO_MAX[int(ld)]

    "======根据地震罕见度，调整地震影响系数最大值======="
    '''
    GB18306-2015 《中国地震动参数区划图》
    基本地震动 basis ground motion 相应于50年超越概率为10%的地震动.
    多遇地震动 frequent ground motion 相应于50年超越概率为63%的地震动。
    罕遇地震动 rare ground motion 相应于50年超越概率为2%的地震动。
    极罕遇地震动 very rare ground motion 相应于50年超越概率为1e-4的地震动。  
    6.2.1 多遇地震动峰值加速度宜按不低于基本地震动峰值加速度1/3倍确定。
    6.2.2 罕遇地震动峰值加速度宜按基本地震动峰值加速度1.6～2.3倍确定。
    6.2.3 极罕遇地震动峰值加速度宜按基本地震动峰值加速度2.7～3.2倍确定。
    地震动峰值加速度 / 重力加速度 = 地震系数  ==========》》》》》》》max(abs(xg))/g = k
    绝对加速度 / 地震动峰值加速度 = 动力放大系数====》》》》 max(abs(xg+x(t))) /  max(abs(xg))  =β（≈2.5，经验）
    地震影响系数定义为 绝对加速度与重力加速度的比值    ==》》》  α = k * β = 2.5 * k 
    因此，地震影响系数为地震动峰值加速度/g的  2.5倍  
    缩放地震动峰值加速度等同于缩放地震影响系数  
    '''
    rarity_dict = {'frequent': 0.33, 'basic': 1, 'rare': 2.3, 'very rare': 3.2}
    arf_max *= rarity_dict[rarity]
    #  根据设计地震分组，场地类别，选择特征周期值
    TG = np.array([[0.2, 0.25, 0.35, 0.45, 0.65],
                   [0.25, 0.3, 0.4, 0.55, 0.75],
                   [0.3, 0.35, 0.45, 0.65, 0.9]])
    if rarity == 'basic' or 'frequent':
        Tg = TG[fz - 1, cd]
    else:  # 计算罕遇地震作用时，特征周期应增加0.05
        Tg = TG[fz - 1, cd] + 0.5

    # 设计谱曲线的参数
    gama = 0.9 + (0.05 - Damp) / (0.3 + 6 * Damp)  # 曲线下降段的衰减指数
    eta1 = 0.02 + (0.05 - Damp) / (4 + 32 * Damp)  # 直线下降段的下降斜率调整系数
    eta2 = 1 + (0.05 - Damp) / (0.08 + 1.6 * Damp)  # 阻尼调整系数
    eta1 = max(0, eta1)
    eta2 = max(eta2, 0.55)
    # 分段位置 T1 T2 T3
    T1 = 0.1
    T2 = Tg  # 特征周期
    T3 = 5 * Tg
    T_jg = T  # 结构周期T

    def f1(T):
        return arf_max * (0.40 + (eta2 - 0.40) / 0.1 * T)

    def f2():
        return eta2 * arf_max

    def f3(T):
        return ((Tg / T) ** gama) * eta2 * arf_max

    def f4(T):
        return (eta2 * 0.2 ** gama - eta1 * (T - 5 * Tg)) * arf_max

    # 第一段 0～T1
    if 0 <= T_jg < T1:
        arf_jg = f1(T_jg)
    # 第二段 T1～T2
    if (T1 <= T_jg) & (T_jg < T2):
        arf_jg = f2()
    # 第三段 T2~T3
    if (T2 <= T_jg) & (T_jg < T3):
        arf_jg = f3(T_jg)
    # 第四段 T3～6.0
    if (T3 <= T_jg) & (T_jg < 6):
        arf_jg = f4(T_jg)
    # 第五段 6.0～
    if 6.0 <= T_jg:
        arf_jg = (eta2 * 0.2 ** gama - eta1 * (6.0 - 5 * Tg)) * arf_max

    # 计算地震影响系数最大值场地调整系数 国标 表5.0.5
    # 6 7 7.5 8 8.25 8.5 >=9
    TZ = np.array([[0.72, 0.74, 0.75, 0.76, 0.805, 0.85, 0.9],
                   [0.8, 0.82, 0.83, 0.85, 0.90, 0.95, 1],
                   [1, 1, 1, 1, 1, 1, 1],
                   [1.3, 1.25, 1.15, 1, 1, 1, 1],
                   [1.25, 1.20, 1.10, 1.00, 0.975, 0.95, 0.9]])
    eta3 = TZ[cd, min(ld, 6)]

    # 反应谱值 拟加速度值
    return arf_jg * eta3  # 返回影响系数  乘上g就是加速度


def designRsFre(T, Damp, ld, cd, fz, rarity):
    """
    # 地震波设计谱函数 用于求取一个频率在设计谱上对应的谱值
    # 横坐标是频率时使用
    # 修改1s以后的曲线 使得反应谱以频率为坐标时更美观
    # 针对电力设施
    :param T: 特征周期序列
    :param Damp: 阻尼比
    :param ld: 地震烈度  6 7 7.5 8 8.25 8.5 9 10
    :param cd: 场地类型 0 1 2 3 4
    :param fz: 场地分组 1 2 3
    :param rarity: 地震罕见度 frequent多遇地震 rare罕遇地震 basic基本地震动 very rare 极罕遇地震动
    :return:
    """
    # 根据抗震设防烈度，选择水平地震影响系数最大值
    arf_jg = 0
    # 水平地震影响系数最大值 6-0.05g 7-0.10g 7.5-0.15g 8-0.20g 8.25-0.25g 8.5-0.30g 9-0.40g 10-0.50g
    ACO_MAX = np.array([0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 1, 1.25])
    arf_max = ACO_MAX[int(ld)]

    "======根据地震罕见度，调整地震影响系数最大值======="
    '''
    GB18306-2015 《中国地震动参数区划图》
    基本地震动 basis ground motion 相应于50年超越概率为10%的地震动.
    多遇地震动 frequent ground motion 相应于50年超越概率为63%的地震动。
    罕遇地震动 rare ground motion 相应于50年超越概率为2%的地震动。
    极罕遇地震动 very rare ground motion 相应于50年超越概率为1e-4的地震动。  
    6.2.1 多遇地震动峰值加速度宜按不低于基本地震动峰值加速度1/3倍确定。
    6.2.2 罕遇地震动峰值加速度宜按基本地震动峰值加速度1.6～2.3倍确定。
    6.2.3 极罕遇地震动峰值加速度宜按基本地震动峰值加速度2.7～3.2倍确定。
    地震动峰值加速度 / 重力加速度 = 地震系数  ==========》》》》》》》max(abs(xg‘’))/g = k
    绝对加速度 / 地震动峰值加速度 = 动力放大系数====》》》》 max(abs(xg''+x(t)'')) /  max(abs(xg‘’))  =β（≈2.5，经验）
    地震影响系数定义为 绝对加速度与重力加速度的比值    ==》》》  α = k * β = 2.5 * k 
    因此，地震影响系数为地震动峰值加速度/g的  2.5倍  
    缩放地震动峰值加速度等同于缩放地震影响系数  
    '''
    rarity_dict = {'frequent': 0.33, 'basic': 1, 'rare': 2.3, 'very rare': 3.2}
    arf_max *= rarity_dict[rarity]
    #  根据设计地震分组，场地类别，选择特征周期值
    TG = np.array([[0.2, 0.25, 0.35, 0.45, 0.65],
                   [0.25, 0.3, 0.4, 0.55, 0.75],
                   [0.3, 0.35, 0.45, 0.65, 0.9]])
    if rarity == 'basic' or 'frequent':
        Tg = TG[fz - 1, cd]
    else:  # 计算罕遇地震作用时，特征周期应增加0.05
        Tg = TG[fz - 1, cd] + 0.5

    # 设计谱曲线的参数
    gama = 0.9 + (0.05 - Damp) / (0.3 + 6 * Damp)  # 曲线下降段的衰减指数
    eta1 = 0.02 + (0.05 - Damp) / (4 + 32 * Damp)  # 直线下降段的下降斜率调整系数
    eta2 = 1 + (0.05 - Damp) / (0.08 + 1.6 * Damp)  # 阻尼调整系数
    eta1 = max(0, eta1)
    eta2 = max(eta2, 0.55)
    # 分段位置 T1 T2 T3
    T1 = 0.1
    T2 = Tg  # 特征周期
    T3 = 5 * Tg
    T_jg = T  # 结构周期T

    def f1(T):
        return arf_max * (0.40 + (eta2 - 0.40) / 0.1 * T)

    def f2():
        return eta2 * arf_max

    def f3(T):
        return ((Tg / T) ** gama) * eta2 * arf_max

    def f4(T):
        return (eta2 * 0.2 ** gama - eta1 * (T - 5 * Tg)) * arf_max

    # 第一段 0～T1
    if 0 <= T_jg < T1:
        arf_jg = f1(T_jg)
    # 第二段 T1～T2
    if (T1 <= T_jg) & (T_jg < T2):
        arf_jg = f2()
    # 第三段 T2~T3
    if (T2 <= T_jg) & (T_jg < T3):
        arf_jg = f3(T_jg)
    # 第四段 T3～6.0
    if (T3 <= T_jg) & (T_jg < 6):
        arf_jg = f3(T_jg)
    # 第五段 6.0～
    if 6.0 <= T_jg:
        arf_jg = f3(T_jg)

    # 计算地震影响系数最大值场地调整系数
    # 6 7 7.5 8 8.25 8.5 >=9
    TZ = np.array([[0.72, 0.74, 0.75, 0.76, 0.805, 0.85, 0.9],
                   [0.8, 0.82, 0.83, 0.85, 0.90, 0.95, 1],
                   [1, 1, 1, 1, 1, 1, 1],
                   [1.3, 1.25, 1.15, 1, 1, 1, 1],
                   [1.25, 1.20, 1.10, 1.00, 0.975, 0.95, 0.9]])
    eta3 = TZ[cd, min(ld, 6)]

    # 反应谱值 拟加速度值
    return arf_jg * eta3  # 返回影响系数  乘上g就是加速度



# # 测试
# if __name__ == '__main__':
#     import matplotlib.pyplot as plt
#     Damp = 0.05
#     ld = 6
#     cd = 2
#     fz = 3
#     rare1 = 'rare'
#     rare2 = 'basic'
#
#     t1 = np.arange(0.1, 3, 0.01)
#     t2 = np.arange(3, 3.6, 0.15)
#     t3 = np.arange(3.6, 5, 0.2)
#     t4 = np.arange(5, 8, 0.25)
#     t5 = np.arange(8, 15, 0.5)
#     t6 = np.arange(15, 34, 1)
#     tt = np.concatenate((t1, t2, t3, t4, t5, t6))
#     # 反应谱频率范围：0.2-34Hz
#     Tyz = np.sort(1. / tt)  # 反应谱横坐标  非均匀
#     b = np.zeros_like(Tyz)
#     b2 = np.zeros_like(Tyz)
#     for i in range(len(Tyz)):
#         b[i] = designRs(Tyz[i], Damp, ld, cd, fz, rare1)
#         b2[i] = designRs(Tyz[i], Damp, ld, cd, fz, rare2)
#
#     plt.subplot(211)
#     plt.plot(Tyz, b)
#     plt.grid()
#     plt.xlabel('T')
#     plt.ylabel('Tyz/g')
#     plt.title('Earthquake influence coefficient for electrics')
#
#     plt.subplot(212)
#     plt.plot(Tyz, b2)
#     plt.grid()
#     plt.xlabel('T')
#     plt.ylabel('Tyz/g')
#     plt.title('Earthquake influence coefficient for buildings')
#     plt.show()
